Loading packages

knitr::opts_chunk$set(echo = TRUE)
pacman::p_load(tidyverse, 
               here,
               metafor,
               emmeans,
               orchaRd, 
               MCMCglmm,
               corpcor)

Loading data and functions

dat <- read_csv(here("Data","Full_data.csv"))
# Load custom function to extract data 
source(here("R/functions.R")) 

Data organisation

Getting effect sizes from function, ‘flipping’ effect sizes so that all effect sizes are higher values = individuals do better and learning/memory, and shifting negative values to possitive as lnRR cannot use negative values

#removing negative values
#dat <- dat %>%
 # mutate(CC_mean = case_when(ES_ID == 36 ~ CC_mean+29.4),
  #       CS_mean = case_when(ES_ID == 36 ~ CC_mean+29.4),
   #      EC_mean = case_when(ES_ID == 36 ~ CC_mean+29.4),
    #     ES_mean = case_when(ES_ID == 36 ~ CC_mean+29.4),
     #    CC_mean = case_when(ES_ID == 37 ~ CC_mean+36.65),
      #   CS_mean = case_when(ES_ID == 37 ~ CC_mean+36.65),
       #  EC_mean = case_when(ES_ID == 37 ~ CC_mean+36.65),
        # ES_mean = case_when(ES_ID == 37 ~ CC_mean+36.65),
         #CC_mean = case_when(ES_ID == 38 ~ CC_mean+80),
         #CS_mean = case_when(ES_ID == 38 ~ CC_mean+80),
         #EC_mean = case_when(ES_ID == 38 ~ CC_mean+80),
         #ES_mean = case_when(ES_ID == 38 ~ CC_mean+80),
         #CC_mean = case_when(ES_ID == 39 ~ CC_mean+33.86),
         #CS_mean = case_when(ES_ID == 39 ~ CC_mean+33.86),
         #EC_mean = case_when(ES_ID == 39 ~ CC_mean+33.86),
         #ES_mean = case_when(ES_ID == 39 ~ CC_mean+33.86))

#dat$CC_mean <- ifelse(dat$ES_ID == 36, dat$CC_mean+29.4, dat$CC_mean)


#Getting effect sizes
effect_size <- effect_set(CC_n = "CC_n", CC_mean = "CC_mean", CC_SD = "CC_SD",
                          EC_n = "EC_n", EC_mean = "EC_mean" , EC_SD ="EC_SD",
                          CS_n = "CS_n", CS_mean = "CS_mean", CS_SD = "CS_SD",
                          ES_n = "ES_n", ES_mean = "ES_mean", ES_SD = "ES_SD",
                          data = dat)
#'pure' effect sizes
effect_size2 <- effect_set2(CC_n = "CC_n", CC_mean = "CC_mean", CC_SD = "CC_SD",
                          EC_n = "EC_n", EC_mean = "EC_mean" , EC_SD ="EC_SD",
                          CS_n = "CS_n", CS_mean = "CS_mean", CS_SD = "CS_SD",
                          ES_n = "ES_n", ES_mean = "ES_mean", ES_SD = "ES_SD",
                          data = dat) 

#Removing missing effect sizes
dim(dat)
full_info <- which(complete.cases(effect_size) == TRUE) 
dat_effect <- cbind(dat, effect_size, effect_size2)
dat <- dat_effect[full_info, ]
names(dat_effect)

# TODO - we need to sort this out a bit more later
#There are no longer any missing effect sizes but is missing info in moderators: do we still need to remove these?
dat <- dat_effect[full_info, ]

dim(dat_effect)
dimentions <- dim(dat) # 7 less

#Flipping 'lower is better' effect sizes
#flipping lnRR for values where higher = worse
dat$lnRR_Ea <- ifelse(dat$Response_direction == 2, dat$lnRR_E*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_E)) # currently NAswhich causes error
dat$lnRR_Sa  <- ifelse(dat$Response_direction == 2, dat$lnRR_S*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_S)) # currently NAswhich causes error
dat$lnRR_ESa <-  ifelse(dat$Response_direction == 2, dat$lnRR_ES*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_ES)) # currently NAswhich causes error

#flipping 'pure effect sizes'
dat$lnRR_E2a <- ifelse(dat$Response_direction == 2, dat$lnRR_E2*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_E2)) # currently NAswhich causes error
dat$lnRR_S2a  <- ifelse(dat$Response_direction == 2, dat$lnRR_S2*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_S2)) # currently NAswhich causes error
dat$lnRR_ES2a <-  ifelse(dat$Response_direction == 2, dat$lnRR_ES2*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_ES2)) # currently NAswhich causes error
dat$lnRR_E3a <-  ifelse(dat$Response_direction == 2, dat$lnRR_E3*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_E3)) # currently NAswhich causes error
dat$lnRR_S3a <-  ifelse(dat$Response_direction == 2, dat$lnRR_S3*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$lnRR_S3)) # currently NAswhich causes error

#flipping SMD
dat$SMD_Ea <- ifelse(dat$Response_direction == 2, dat$SMD_E*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$SMD_E)) # currently NAswhich causes error
dat$SMD_Sa  <- ifelse(dat$Response_direction == 2, dat$SMD_S*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$SMD_S)) # currently NAswhich causes error
dat$SMD_ESa <-  ifelse(dat$Response_direction == 2, dat$SMD_ES*-1,ifelse(is.na(dat$Response_direction) == TRUE, NA, dat$SMD_ES))

#rounding down integer sample sizes 
dat$CC_n <- floor(dat$CC_n)
dat$EC_n <- floor(dat$EC_n)
dat$CS_n <- floor(dat$CS_n)
dat$ES_n <- floor(dat$CS_n)

Assigning human readable words to moderators

dat <- dat %>% mutate(Type_learning = case_when(Type_learning == 1 ~ "Habituation",
                                                Type_learning == 2 ~ "Conditioning",
                                                Type_learning == 3 ~ "Recognition", 
                                                Type_learning == 4 ~ "Unclear"),
                      Learning_vs_memory = case_when(Learning_vs_memory == 1 ~ "Learning",
                                                     Learning_vs_memory == 2 ~ "Memory", 
                                                     Learning_vs_memory == 1 ~ "Unclear"),
                      Appetitive_vs_aversive = case_when(Appetitive_vs_aversive == 1 ~"Appetitive",
                                                         Appetitive_vs_aversive == 2 ~ "Aversive",
                                                         Appetitive_vs_aversive == 3 ~ "Not applicable",
                                                         Appetitive_vs_aversive == 4 ~ "Unclear"),
                      Type_stress_exposure = case_when(Type_stress_exposure == 1 ~ "Density",
                                                       Type_stress_exposure == 2 ~ "Scent",
                                                       Type_stress_exposure == 3 ~ "Shock",
                                                       Type_stress_exposure == 4 ~ "Exertion",
                                                       Type_stress_exposure == 5 ~ "Restraint",
                                                       Type_stress_exposure == 6 ~ "MS",
                                                       Type_stress_exposure == 7 ~ "Circadian rhythm",
                                                       Type_stress_exposure == 8 ~ "Noise",
                                                       Type_stress_exposure == 9 ~ "Other",
                                                       Type_stress_exposure == 10 ~ "Combination",
                                                       Type_stress_exposure == 11 ~ "unclear"), 
                      Age_stress_exposure = case_when(Age_stress_exposure == 1 ~ "Prenatal",
                                                      Age_stress_exposure == 2 ~ "Juvenile",
                                                      Age_stress_exposure == 3 ~ "Adult",
                                                      Age_stress_exposure == 4 ~ "Unclear"),
                      Stress_duration = case_when(Stress_duration == 1 ~ "Acute",
                                                  Stress_duration == 2 ~ "Chronic",
                                                  Stress_duration == 3 ~ "Intermittent",
                                                  Stress_duration == 4 ~ "Unclear"),
                      EE_social = case_when(EE_social == 1 ~ "Social",
                                            EE_social== 2 ~ "Non-social", 
                                            EE_social == 3 ~ "Unclear"), 
                      EE_exercise = case_when(EE_exercise == 1 ~ "Exercise", 
                                              EE_exercise == 2 ~ "No exercise"),
                      Age_EE_exposure = case_when(Age_EE_exposure == 1 ~ "Prenatal", 
                                                  Age_EE_exposure == 2 ~ "Juvenile",
                                                  Age_EE_exposure == 3 ~ "Adult", 
                                                  Age_EE_exposure == 4 ~ "Unclear"))

Making Variance VCV

Used to control for trait relatedness within individuals

##adding data to empty VCV matrix
VCV_E <- make_VCV_matrix(data = dat, V = "lnRRV_E", cluster = "Study_ID", obs = "ES_ID")
VCV_S <- make_VCV_matrix(data = dat, V = "lnRRV_S", cluster = "Study_ID", obs = "ES_ID")
VCV_ES <- make_VCV_matrix(data = dat, V = "lnRRV_ES", cluster = "Study_ID", obs = "ES_ID")

is.positive.definite(VCV_E)
isSymmetric(VCV_E)
is.positive.definite(VCV_S)
isSymmetric(VCV_S)
is.positive.definite(VCV_ES)
isSymmetric(VCV_ES)

VCV_Ea <- make_VCV_matrix(data = dat, V = "SMDV_E", cluster = "Study_ID", obs = "ES_ID")
VCV_Sa <- make_VCV_matrix(data = dat, V = "SMDV_S", cluster = "Study_ID", obs = "ES_ID")
VCV_ESa <- make_VCV_matrix(data = dat, V = "SMDV_ES", cluster = "Study_ID", obs = "ES_ID")

is.positive.definite(VCV_Ea)
isSymmetric(VCV_Ea)
is.positive.definite(VCV_Sa)
isSymmetric(VCV_Sa)
is.positive.definite(VCV_ESa)
isSymmetric(VCV_ESa)

Examining the group of large/positive effect sizes

Note that I checked the top four most positive interaction effect sizes and the data looks correct. Highest ES is study 27, effect size 88, where Maternally separated EE individuals freeze significantly more when given a cue from inhibitory avoidance. Im not sure why the second highest one (ES 16) is so high as the control indivs do the best. EE indivs do best in the 3rd highest ES (ES 10) so not sure about this one either? The fourth highest (ES 12) seems correct

largeES_S <- dat %>%
  slice_max(lnRR_S, n = 10)

largeES_E <- dat %>%
  slice_max(lnRR_E, n = 10)

largeES_ES <-dat %>%
  slice_max(lnRR_ES, n = 10)

#most of the largest effect sizes come from different papers

Modelling with lnRR

Things todo: - Incorporate strain as random effect - Consider including VCV - Check outlier - Remove moderator levels with k < 5

Stress

Intercept model

Learning and memory are significantly reduced due to stress. High heterogeneity

mod_S0 <- rma.mv(yi = lnRR_Sa, V = VCV_S, random = list(~1|Study_ID,
                                                          ~1|ES_ID,
                                                          ~1|Strain),
                 test = "t", 
                 data = dat)

summary(mod_S0) 
## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -10.5945   21.1890   29.1890   39.2325   29.6542   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0080  0.0892     30     no  Study_ID 
## sigma^2.2  0.0358  0.1891     92     no     ES_ID 
## sigma^2.3  0.0001  0.0086      6     no    Strain 
## 
## Test for Heterogeneity:
## Q(df = 91) = 863.4222, p-val < .0001
## 
## Model Results:
## 
## estimate      se     tval  df    pval    ci.lb    ci.ub 
##  -0.0903  0.0331  -2.7240  91  0.0077  -0.1561  -0.0244  ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_S0) 
##    I2_total I2_Study_ID    I2_ES_ID   I2_Strain 
## 0.928503647 0.168668323 0.758268848 0.001566476
funnel(mod_S0)

orchard_plot(mod_S0, mod = "Int", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Publication Bias

… I was going to run this using the same code as last meta-analysis (MCMCglmm) but thought there might be an easier way?

#need to remove NAs

dat_na_S <- dat %>%
  drop_na(c(lnRR_Sa, lnRRV_S))

dat_na_E <- dat  %>%
  drop_na(c(lnRR_Ea, lnRRV_E))

dat_na_ES <- dat  %>%
  drop_na(c(lnRR_ESa, lnRRV_ES))

#Stress
#TODO check with Shinichi if there is an easier way to do PB (eggers regression, funnel asymmetry)

Meta-regression

Type of learning

The type of learning/memory response

dat$Type_learning<-as.factor(dat$Type_learning)

dat1 <- filter(dat, Type_learning %in% c("Recognition", "Habituation", "Conditioning"))
VCV_S1 <- make_VCV_matrix(data = dat1, V = "lnRRV_S", cluster = "Study_ID", obs = "ES_ID")


mod_S1 <- rma.mv(yi = lnRR_Sa, V = VCV_S1, mod = ~Type_learning-1, random =   list(~1|Study_ID,
                                                                                    ~1|ES_ID,
                                                                                    ~1|Strain),
                 test = "t",
                 data = dat1)

summary(mod_S1)
## 
## Multivariate Meta-Analysis Model (k = 91; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -4.8116    9.6232   21.6232   36.4872   22.6602   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0196  0.1400     30     no  Study_ID 
## sigma^2.2  0.0192  0.1385     91     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 88) = 498.8005, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 88) = 5.2846, p-val = 0.0021
## 
## Model Results:
## 
##                            estimate      se     tval  df    pval    ci.lb 
## Type_learningConditioning   -0.1110  0.0387  -2.8663  88  0.0052  -0.1880 
## Type_learningHabituation    -0.2548  0.0816  -3.1212  88  0.0024  -0.4170 
## Type_learningRecognition     0.0096  0.0571   0.1686  88  0.8665  -0.1038 
##                              ci.ub 
## Type_learningConditioning  -0.0340  ** 
## Type_learningHabituation   -0.0926  ** 
## Type_learningRecognition    0.1230     
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S1) 
##   R2_marginal R2_coditional 
##    0.08953573    1.00000000
# Orchard plot 
orchard_plot(mod_S1, mod = "Type_learning", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Learning vs Memory

Is the assay broadly measuring learning or memory?

mod_S2 <-  rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID,
                                                                                        ~1|ES_ID,
                                                                                        ~1|Strain),
                  test = "t",
                  data = dat)

summary(mod_S2) 
## 
## Multivariate Meta-Analysis Model (k = 85; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -10.8490   21.6980   31.6980   43.7922   32.4772   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0184  0.1357     30     no  Study_ID 
## sigma^2.2  0.0279  0.1669     85     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 83) = 623.2963, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 83) = 2.7188, p-val = 0.0718
## 
## Model Results:
## 
##                             estimate      se     tval  df    pval    ci.lb 
## Learning_vs_memoryLearning   -0.0342  0.0509  -0.6715  83  0.5037  -0.1354 
## Learning_vs_memoryMemory     -0.0922  0.0400  -2.3073  83  0.0235  -0.1718 
##                               ci.ub 
## Learning_vs_memoryLearning   0.0670    
## Learning_vs_memoryMemory    -0.0127  * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S2) 
##   R2_marginal R2_coditional 
##    0.01603098    1.00000000
# Orchard plot 
orchard_plot(mod_S2, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Appetitive vs aversive

The type of cue used

dat2 <- filter(dat, Appetitive_vs_aversive %in% c("Appetitive", "Aversive", "Not applicable"))
VCV_S2 <- make_VCV_matrix(data = dat2, V = "lnRRV_S", cluster = "Study_ID", obs = "ES_ID")

mod_S3 <- rma.mv(yi = lnRR_Sa, V = VCV_S2, mod = ~ Appetitive_vs_aversive-1, random = list(~1|Study_ID,
                                                                                            ~1|ES_ID,
                                                                                            ~1|Strain),
                 test = "t",
                 data = dat2)

summary(mod_S3)
## 
## Multivariate Meta-Analysis Model (k = 91; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -6.5287   13.0575   25.0575   39.9215   26.0945   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0145  0.1203     30     no  Study_ID 
## sigma^2.2  0.0256  0.1601     91     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 88) = 506.6327, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 88) = 3.8651, p-val = 0.0120
## 
## Model Results:
## 
##                                       estimate      se     tval  df    pval 
## Appetitive_vs_aversiveAppetitive       -0.1781  0.0756  -2.3563  88  0.0207 
## Appetitive_vs_aversiveAversive         -0.1065  0.0429  -2.4801  88  0.0150 
## Appetitive_vs_aversiveNot applicable   -0.0320  0.0508  -0.6306  88  0.5300 
##                                         ci.lb    ci.ub 
## Appetitive_vs_aversiveAppetitive      -0.3283  -0.0279  * 
## Appetitive_vs_aversiveAversive        -0.1919  -0.0212  * 
## Appetitive_vs_aversiveNot applicable  -0.1329   0.0689    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S3) 
##   R2_marginal R2_coditional 
##    0.05569185    1.00000000
# Orchard plot 
orchard_plot(mod_S3, mod = "Appetitive_vs_aversive", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Type of stress

The type of stress manipulation

dat3 <- filter(dat, Type_stress_exposure %in% c("Restraint", "Noise", "MS", "Combination"))
VCV_S3 <- make_VCV_matrix(data = dat3, V = "lnRRV_S", cluster = "Study_ID", obs = "ES_ID")

mod_S4 <- rma.mv(yi = lnRR_Sa, V = VCV_S3, mod = ~Type_stress_exposure-1, random = list(~1|Study_ID,
                                                                                         ~1|ES_ID,
                                                                                         ~1|Strain),
                 test = "t",
                 data = dat3)
summary(mod_S4) 
## 
## Multivariate Meta-Analysis Model (k = 87; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -9.2101   18.4202   32.4202   49.3521   33.9135   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0080  0.0896     26     no  Study_ID 
## sigma^2.2  0.0373  0.1932     87     no     ES_ID 
## sigma^2.3  0.0000  0.0000      4     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 83) = 764.1363, p-val < .0001
## 
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 83) = 2.7813, p-val = 0.0319
## 
## Model Results:
## 
##                                  estimate      se     tval  df    pval    ci.lb 
## Type_stress_exposureCombination   -0.0809  0.0851  -0.9511  83  0.3443  -0.2502 
## Type_stress_exposureMS            -0.0507  0.0508  -0.9983  83  0.3210  -0.1518 
## Type_stress_exposureNoise         -0.0739  0.0905  -0.8164  83  0.4166  -0.2539 
## Type_stress_exposureRestraint     -0.2026  0.0692  -2.9253  83  0.0044  -0.3403 
##                                    ci.ub 
## Type_stress_exposureCombination   0.0883     
## Type_stress_exposureMS            0.0503     
## Type_stress_exposureNoise         0.1061     
## Type_stress_exposureRestraint    -0.0648  ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S4)
##   R2_marginal R2_coditional 
##     0.0690348     1.0000000
# Orchard plot 
orchard_plot(mod_S4, mod = "Type_stress_exposure", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Age of stress

The age at which the individuals were exposed to the stressor. Note there are a lot of ‘unkown’ age as authors only report PND which needs to be researched. I’m wondering if this also needs an ‘adolescence’ category as this seesm to be popular in rodent research

#TODO these ages need to be recategoried and analysis run again

mod_S5 <-rma.mv(yi = lnRR_Sa, V = VCV_S, mod = ~Age_stress_exposure-1, random = list(~1|Study_ID,
                                                                                       ~1|ES_ID,
                                                                                       ~1|Strain),
                test = "t",
                data = dat)
summary(mod_S5) 
## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -9.2914   18.5827   32.5827   49.9241   33.9827   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0014  0.0372     30     no  Study_ID 
## sigma^2.2  0.0381  0.1953     92     no     ES_ID 
## sigma^2.3  0.0032  0.0565      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 88) = 800.8761, p-val < .0001
## 
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 88) = 1.9490, p-val = 0.1094
## 
## Model Results:
## 
##                              estimate      se     tval  df    pval    ci.lb 
## Age_stress_exposureAdult      -0.1467  0.0857  -1.7123  88  0.0904  -0.3171 
## Age_stress_exposureJuvenile   -0.0012  0.0543  -0.0214  88  0.9829  -0.1090 
## Age_stress_exposurePrenatal   -0.1276  0.0873  -1.4617  88  0.1474  -0.3011 
## Age_stress_exposureUnclear    -0.1588  0.0874  -1.8172  88  0.0726  -0.3325 
##                               ci.ub 
## Age_stress_exposureAdult     0.0236  . 
## Age_stress_exposureJuvenile  0.1067    
## Age_stress_exposurePrenatal  0.0459    
## Age_stress_exposureUnclear   0.0149  . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S5) 
##   R2_marginal R2_coditional 
##     0.1009404     0.9327030
# Orchard plot 
orchard_plot(mod_S5, mod = "Age_stress_exposure", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Acute vs chronic stress

How long was the stress applied for (chronic = every day for 7 days or more)? This has the highest marginal R2

dat4 <- filter(dat, Stress_duration %in% c("Chronic", "Acute"))
VCV_S4 <- make_VCV_matrix(data = dat3, V = "lnRRV_S", cluster = "Study_ID", obs = "ES_ID")

mod_S6 <-rma.mv(yi = lnRR_Sa, V = VCV_S4, mod = ~Stress_duration-1, random = list(~1|Study_ID,
                                                                                   ~1|ES_ID,
                                                                                   ~1|Strain),
                test = "t",
                data = dat4)
summary(mod_S6) 
## 
## Multivariate Meta-Analysis Model (k = 87; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -20.6635   41.3269   51.3269   63.5402   52.0864   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0085  0.0922     27     no  Study_ID 
## sigma^2.2  0.0700  0.2645     87     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 85) = 2809.8421, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 85) = 6.9113, p-val = 0.0017
## 
## Model Results:
## 
##                         estimate      se     tval  df    pval    ci.lb    ci.ub 
## Stress_durationAcute      0.1099  0.0873   1.2581  85  0.2118  -0.0638   0.2835 
## Stress_durationChronic   -0.1569  0.0453  -3.4622  85  0.0008  -0.2470  -0.0668 
##  
## Stress_durationAcute 
## Stress_durationChronic  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_S6) 
##   R2_marginal R2_coditional 
##      0.135429      1.000000
# Orchard plot 
orchard_plot(mod_S6, mod = "Stress_duration", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Environmental enrichment

Intercept model

Learning and memory are significantly improved when housed with environmnetal enrichment

mod_E0 <- rma.mv(yi = lnRR_Ea, V = VCV_E, random = list(~1|Study_ID,
                                                          ~1|ES_ID,
                                                          ~1|Strain),
                 test = "t",
                 data = dat)
summary(mod_E0) 
## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -13.9611   27.9223   35.9223   45.9657   36.3874   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0049  0.0698     30     no  Study_ID 
## sigma^2.2  0.0428  0.2070     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Heterogeneity:
## Q(df = 91) = 914.4170, p-val < .0001
## 
## Model Results:
## 
## estimate      se    tval  df    pval   ci.lb   ci.ub 
##   0.1891  0.0321  5.8832  91  <.0001  0.1253  0.2530  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_E0) 
##     I2_total  I2_Study_ID     I2_ES_ID    I2_Strain 
## 9.339783e-01 9.534725e-02 8.386310e-01 1.741556e-09
funnel(mod_E0)

#trying orchard plot

orchard_plot(mod_E0, mod = "Int", xlab = "lnRR", alpha=0.4) +  # Orchard plot 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5)+ # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2)+ # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_colour_manual(values = "darkorange")+ # change colours
  scale_fill_manual(values="darkorange")+ 
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Meta-regression

Type of learning

The type of learning/memory response

VCV_E1 <- make_VCV_matrix(data = dat1, V = "lnRRV_E", cluster = "Study_ID", obs = "ES_ID")

mod_E1 <- rma.mv(yi = lnRR_Ea, V = VCV_E1, mod = ~Type_learning-1, random = list(~1|Study_ID,
                                                                                  ~1|ES_ID,
                                                                                  ~1|Strain),
                 test = "t",
                 data = dat1)

summary(mod_E1)
## 
## Multivariate Meta-Analysis Model (k = 91; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -11.7462   23.4925   35.4925   50.3565   36.5295   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0145  0.1205     30     no  Study_ID 
## sigma^2.2  0.0346  0.1859     91     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 88) = 823.6996, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 88) = 12.2253, p-val < .0001
## 
## Model Results:
## 
##                            estimate      se    tval  df    pval    ci.lb 
## Type_learningConditioning    0.2415  0.0403  5.9907  88  <.0001   0.1614 
## Type_learningHabituation     0.1760  0.0948  1.8565  88  0.0667  -0.0124 
## Type_learningRecognition     0.0534  0.0654  0.8167  88  0.4163  -0.0765 
##                             ci.ub 
## Type_learningConditioning  0.3217  *** 
## Type_learningHabituation   0.3643    . 
## Type_learningRecognition   0.1833      
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E1) 
##   R2_marginal R2_coditional 
##    0.08695344    1.00000000
# Orchard plot 
orchard_plot(mod_E1, mod = "Type_learning", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Learning vs Memory

Is the assay broadly measuring learning or memory?

mod_E2 <-  rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID,
                                                                                        ~1|ES_ID,
                                                                                        ~1|Strain),
                  test = "t",
                  data = dat)

summary(mod_E2) 
## 
## Multivariate Meta-Analysis Model (k = 85; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -8.9742   17.9484   27.9484   40.0426   28.7276   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0074  0.0860     30     no  Study_ID 
## sigma^2.2  0.0291  0.1707     85     no     ES_ID 
## sigma^2.3  0.0006  0.0255      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 83) = 580.2953, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 83) = 13.4308, p-val < .0001
## 
## Model Results:
## 
##                             estimate      se    tval  df    pval   ci.lb 
## Learning_vs_memoryLearning    0.1991  0.0505  3.9454  83  0.0002  0.0987 
## Learning_vs_memoryMemory      0.1837  0.0392  4.6870  83  <.0001  0.1057 
##                              ci.ub 
## Learning_vs_memoryLearning  0.2994  *** 
## Learning_vs_memoryMemory    0.2616  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E2) 
##   R2_marginal R2_coditional 
##   0.001428381   0.982584068
# Orchard plot 
orchard_plot(mod_E2, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Appetitive vs aversive

The type of cue used

VCV_E2 <- make_VCV_matrix(data = dat2, V = "lnRRV_E", cluster = "Study_ID", obs = "ES_ID")
mod_E3 <- rma.mv(yi = lnRR_Ea, V = VCV_E2, mod = ~ Appetitive_vs_aversive-1, random = list(~1|Study_ID,
                                                                                            ~1|ES_ID,
                                                                                            ~1|Strain),
                 test = "t",
                 data = dat2)

summary(mod_E3)
## 
## Multivariate Meta-Analysis Model (k = 91; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
##  -9.0470   18.0940   30.0940   44.9581   31.1311   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0132  0.1148     30     no  Study_ID 
## sigma^2.2  0.0308  0.1755     91     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 88) = 638.2083, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 88) = 14.8231, p-val < .0001
## 
## Model Results:
## 
##                                       estimate      se    tval  df    pval 
## Appetitive_vs_aversiveAppetitive        0.1895  0.0762  2.4855  88  0.0148 
## Appetitive_vs_aversiveAversive          0.2719  0.0440  6.1808  88  <.0001 
## Appetitive_vs_aversiveNot applicable    0.0623  0.0528  1.1792  88  0.2415 
##                                         ci.lb   ci.ub 
## Appetitive_vs_aversiveAppetitive       0.0380  0.3410    * 
## Appetitive_vs_aversiveAversive         0.1845  0.3594  *** 
## Appetitive_vs_aversiveNot applicable  -0.0427  0.1673      
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E3) 
##   R2_marginal R2_coditional 
##     0.1467484     1.0000000
# Orchard plot 
orchard_plot(mod_E3, mod = "Appetitive_vs_aversive", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Social enrichment

Does EE also include a manipulation of social environment? Note that we excluded any studies that exclusively used social enrichment.s

dat5 <- filter(dat, EE_social %in% c("Social", "Non-social"))
VCV_E5 <- make_VCV_matrix(data = dat5, V = "lnRRV_E", cluster = "Study_ID", obs = "ES_ID")
  
mod_E4<- rma.mv(yi = lnRR_Ea, V = VCV_E5, mod = ~EE_social-1, random = list(~1|Study_ID, 
                                                                             ~1|ES_ID,
                                                                             ~1|Strain),
                test = "t",data = dat5)

summary(mod_E4)
## 
## Multivariate Meta-Analysis Model (k = 90; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -13.6466   27.2932   37.2932   49.6798   38.0249   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0056  0.0747     29     no  Study_ID 
## sigma^2.2  0.0432  0.2078     90     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 88) = 861.6301, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 88) = 16.1474, p-val < .0001
## 
## Model Results:
## 
##                      estimate      se    tval  df    pval   ci.lb   ci.ub 
## EE_socialNon-social    0.1360  0.0535  2.5402  88  0.0128  0.0296  0.2423    * 
## EE_socialSocial        0.2148  0.0423  5.0835  88  <.0001  0.1308  0.2988  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E4) 
##   R2_marginal R2_coditional 
##    0.02972383    1.00000000
# Orchard plot 
orchard_plot(mod_E4, mod = "EE_social", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Exercise enrichment

Does the form of enrichment include exercise through a running wheel or treadmill?

mod_E5<- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~EE_exercise-1, random = list(~1|Study_ID,
                                                                               ~1|ES_ID,
                                                                               ~1|Strain),
                test = "t",
                data = dat)

summary(mod_E5)
## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -14.1841   28.3682   38.3682   50.8673   39.0825   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0055  0.0744     30     no  Study_ID 
## sigma^2.2  0.0431  0.2077     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 90) = 889.4288, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 16.8895, p-val < .0001
## 
## Model Results:
## 
##                         estimate      se    tval  df    pval   ci.lb   ci.ub 
## EE_exerciseExercise       0.1868  0.0396  4.7137  90  <.0001  0.1081  0.2656 
## EE_exerciseNo exercise    0.1953  0.0574  3.3999  90  0.0010  0.0812  0.3094 
##  
## EE_exerciseExercise     *** 
## EE_exerciseNo exercise   ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E5) 
##   R2_marginal R2_coditional 
##  0.0003377362  0.9999999972
# Orchard plot 
orchard_plot(mod_E5, mod = "EE_exercise", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Age of enrichment

The age at which the individuals were exposed to environmental enrichment. This needs to be redone with new age categories

mod_E6 <- rma.mv(yi = lnRR_Ea, V = VCV_E, mod = ~Age_EE_exposure-1, random = list(~1|Study_ID,
                                                                                    ~1|ES_ID,
                                                                                    ~1|Strain),
                 test = "t",
                 data = dat)

summary(mod_E6) 
## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -12.7731   25.5462   37.5462   52.4781   38.5706   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0049  0.0699     30     no  Study_ID 
## sigma^2.2  0.0437  0.2090     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 89) = 882.6717, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 12.5329, p-val < .0001
## 
## Model Results:
## 
##                          estimate      se    tval  df    pval    ci.lb   ci.ub 
## Age_EE_exposureAdult       0.1838  0.0648  2.8351  89  0.0057   0.0550  0.3126 
## Age_EE_exposureJuvenile    0.0082  0.1072  0.0767  89  0.9390  -0.2047  0.2212 
## Age_EE_exposureUnclear     0.2164  0.0398  5.4365  89  <.0001   0.1373  0.2955 
##  
## Age_EE_exposureAdult      ** 
## Age_EE_exposureJuvenile 
## Age_EE_exposureUnclear   *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_E6) 
##   R2_marginal R2_coditional 
##    0.07138796    1.00000000
# Orchard plot 
orchard_plot(mod_E6, mod = "Age_EE_exposure", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Interaction of stress and EE

Intercept

Enriched and stress animals are better at learning and memory. TODO: It looks like there is a large but low precision outlier. Should potentially remove?

mod_ES0 <- rma.mv(yi = lnRR_ESa, V = VCV_ES, random = list(~1|Study_ID,
                                                             ~1|ES_ID,
                                                             ~1|Strain),
                  test = "t", 
                  data = dat)

summary(mod_ES0) 
## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -45.8577   91.7153   99.7153  109.7588  100.1804   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0347  0.1862     30     no  Study_ID 
## sigma^2.2  0.0269  0.1640     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Heterogeneity:
## Q(df = 91) = 322.9690, p-val < .0001
## 
## Model Results:
## 
## estimate      se    tval  df    pval   ci.lb   ci.ub 
##   0.1288  0.0496  2.5996  91  0.0109  0.0304  0.2273  * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_ES0) 
##     I2_total  I2_Study_ID     I2_ES_ID    I2_Strain 
## 8.202982e-01 4.617725e-01 3.585257e-01 1.352473e-08
funnel(mod_ES0)

orchard_plot(mod_ES0, mod = "Int", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

#running the model again with the extreme outlier removed
dat_outliers <- dat %>%
  filter(ES_ID != 88)

VCV_ES_outliers <- make_VCV_matrix(data = dat_outliers, V = "lnRRV_ES", cluster = "Study_ID", obs = "ES_ID")

mod_ESoutlier <- rma.mv(yi = lnRR_ESa, V = VCV_ES_outliers, random = list(~1|Study_ID,
                                                             ~1|ES_ID,
                                                             ~1|Strain),
                  test = "t", 
                  data = dat_outliers)

summary(mod_ES0)
## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -45.8577   91.7153   99.7153  109.7588  100.1804   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0347  0.1862     30     no  Study_ID 
## sigma^2.2  0.0269  0.1640     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Heterogeneity:
## Q(df = 91) = 322.9690, p-val < .0001
## 
## Model Results:
## 
## estimate      se    tval  df    pval   ci.lb   ci.ub 
##   0.1288  0.0496  2.5996  91  0.0109  0.0304  0.2273  * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
orchard_plot(mod_ESoutlier, mod = "Int", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13))   

Meta-regression

Type of learning

The type of learning/memory response

VCV_ES1 <- make_VCV_matrix(data = dat1, V = "lnRRV_ES", cluster = "Study_ID", obs = "ES_ID")

mod_ES1 <- rma.mv(yi = lnRR_ESa, V = VCV_ES1, mod = ~Type_learning-1, random = list(~1|Study_ID,
                                                                                    ~1|ES_ID,
                                                                                    ~1|Strain),
                  test = "t",
                  data = dat1)

summary(mod_ES1)
## 
## Multivariate Meta-Analysis Model (k = 91; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -42.8899   85.7798   97.7798  112.6438   98.8169   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0390  0.1974     30     no  Study_ID 
## sigma^2.2  0.0194  0.1392     91     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 88) = 306.3683, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 88) = 5.3955, p-val = 0.0019
## 
## Model Results:
## 
##                            estimate      se     tval  df    pval    ci.lb 
## Type_learningConditioning    0.1874  0.0550   3.4087  88  0.0010   0.0782 
## Type_learningHabituation     0.1787  0.1069   1.6724  88  0.0980  -0.0336 
## Type_learningRecognition    -0.0429  0.0740  -0.5800  88  0.5634  -0.1900 
##                             ci.ub 
## Type_learningConditioning  0.2967  *** 
## Type_learningHabituation   0.3911    . 
## Type_learningRecognition   0.1041      
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES1) 
##   R2_marginal R2_coditional 
##     0.1063433     1.0000000
# Orchard plot 
orchard_plot(mod_ES1, mod = "Type_learning", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Learning vs Memory

Is the assay broadly measuring learning or memory?

mod_ES2 <-  rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Learning_vs_memory-1, random = list(~1|Study_ID, 
                                                                                           ~1|ES_ID,
                                                                                           ~1|Strain),
                   test = "t",
                   data = dat)

summary(mod_ES2) 
## 
## Multivariate Meta-Analysis Model (k = 85; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -47.7835   95.5669  105.5669  117.6611  106.3461   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0345  0.1857     30     no  Study_ID 
## sigma^2.2  0.0313  0.1769     85     no     ES_ID 
## sigma^2.3  0.0000  0.0001      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 83) = 313.8032, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 83) = 3.7300, p-val = 0.0281
## 
## Model Results:
## 
##                             estimate      se    tval  df    pval    ci.lb 
## Learning_vs_memoryLearning    0.1817  0.0694  2.6157  83  0.0106   0.0435 
## Learning_vs_memoryMemory      0.1066  0.0540  1.9724  83  0.0519  -0.0009 
##                              ci.ub 
## Learning_vs_memoryLearning  0.3198  * 
## Learning_vs_memoryMemory    0.2141  . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES2) 
##   R2_marginal R2_coditional 
##     0.0187792     0.9999999
# Orchard plot 
orchard_plot(mod_ES2, mod = "Learning_vs_memory", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Appetitive vs aversive

The type of cue used

VCV_ES2 <- make_VCV_matrix(data = dat2, V = "lnRRV_ES", cluster = "Study_ID", obs = "ES_ID")
mod_ES3 <- rma.mv(yi = lnRR_ESa, V = VCV_ES2, mod = ~ Appetitive_vs_aversive-1, random = list(~1|Study_ID,
                                                                                               ~1|ES_ID,
                                                                                               ~1|Strain),
                  test = "t",
                  data = dat2)

summary(mod_ES3)
## 
## Multivariate Meta-Analysis Model (k = 91; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -43.0121   86.0243   98.0243  112.8883   99.0613   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0327  0.1809     30     no  Study_ID 
## sigma^2.2  0.0218  0.1476     91     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 88) = 281.9381, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 88) = 5.0524, p-val = 0.0028
## 
## Model Results:
## 
##                                       estimate      se    tval  df    pval 
## Appetitive_vs_aversiveAppetitive        0.1612  0.1081  1.4909  88  0.1396 
## Appetitive_vs_aversiveAversive          0.1986  0.0594  3.3409  88  0.0012 
## Appetitive_vs_aversiveNot applicable    0.0126  0.0652  0.1929  88  0.8475 
##                                         ci.lb   ci.ub 
## Appetitive_vs_aversiveAppetitive      -0.0537  0.3761     
## Appetitive_vs_aversiveAversive         0.0805  0.3168  ** 
## Appetitive_vs_aversiveNot applicable  -0.1169  0.1421     
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES3) 
##   R2_marginal R2_coditional 
##    0.09870464    1.00000000
# Orchard plot 
orchard_plot(mod_ES3, mod = "Appetitive_vs_aversive", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Type of stress

The type of stress manipulation

VCV_ES3 <- make_VCV_matrix(data = dat3, V = "lnRRV_ES", cluster = "Study_ID", obs = "ES_ID")
mod_ES4 <- rma.mv(yi = lnRR_ESa, V = VCV_ES3, mod = ~Type_stress_exposure-1, random = list(~1|Study_ID,
                                                                                            ~1|ES_ID,
                                                                                            ~1|Strain),
                  test = "t",
                  data = dat3)
summary(mod_ES4)
## 
## Multivariate Meta-Analysis Model (k = 87; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -41.9621   83.9243   97.9243  114.8562   99.4176   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0469  0.2167     26     no  Study_ID 
## sigma^2.2  0.0274  0.1655     87     no     ES_ID 
## sigma^2.3  0.0000  0.0000      4     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 83) = 292.2191, p-val < .0001
## 
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 83) = 1.3135, p-val = 0.2718
## 
## Model Results:
## 
##                                  estimate      se    tval  df    pval    ci.lb 
## Type_stress_exposureCombination    0.0995  0.1398  0.7121  83  0.4784  -0.1785 
## Type_stress_exposureMS             0.1167  0.0853  1.3684  83  0.1749  -0.0529 
## Type_stress_exposureNoise          0.1709  0.1592  1.0735  83  0.2862  -0.1458 
## Type_stress_exposureRestraint      0.1480  0.1128  1.3122  83  0.1931  -0.0764 
##                                   ci.ub 
## Type_stress_exposureCombination  0.3775    
## Type_stress_exposureMS           0.2863    
## Type_stress_exposureNoise        0.4876    
## Type_stress_exposureRestraint    0.3724    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES4)
##   R2_marginal R2_coditional 
##   0.007436653   0.999999995
# Orchard plot 
orchard_plot(mod_ES4, mod = "Type_stress_exposure", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Age of stress

The age at which the individuals were exposed to the stressor. Note there are a lot of ‘unkown’ age as authors only report PND which needs to be researched. I’m wondering if this also needs an ‘adolescence’ category as this seesm to be popular in rodent research This nees to be redone after age categories are re done

mod_ES5 <-rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~Age_stress_exposure-1, random = list(~1|Study_ID,
                                                                                          ~1|ES_ID,
                                                                                          ~1|Strain),
                 test = "t",
                 data = dat)
summary(mod_ES5) 
## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -42.6838   85.3676   99.3676  116.7089  100.7676   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0238  0.1542     30     no  Study_ID 
## sigma^2.2  0.0259  0.1610     92     no     ES_ID 
## sigma^2.3  0.0131  0.1146      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 88) = 282.2171, p-val < .0001
## 
## Test of Moderators (coefficients 1:4):
## F(df1 = 4, df2 = 88) = 2.2788, p-val = 0.0671
## 
## Model Results:
## 
##                              estimate      se    tval  df    pval    ci.lb 
## Age_stress_exposureAdult       0.0739  0.1358  0.5438  88  0.5880  -0.1960 
## Age_stress_exposureJuvenile    0.0472  0.0941  0.5008  88  0.6177  -0.1399 
## Age_stress_exposurePrenatal    0.3611  0.1492  2.4204  88  0.0176   0.0646 
## Age_stress_exposureUnclear     0.2506  0.1345  1.8624  88  0.0659  -0.0168 
##                               ci.ub 
## Age_stress_exposureAdult     0.3437    
## Age_stress_exposureJuvenile  0.2342    
## Age_stress_exposurePrenatal  0.6576  * 
## Age_stress_exposureUnclear   0.5179  . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES5) 
##   R2_marginal R2_coditional 
##     0.1611826     0.8245990
# Orchard plot 
orchard_plot(mod_ES5, mod = "Age_stress_exposure", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Acute vs chronic stress

How long was the stress applied for (chronic = every day for 7 days or more)? This has the highest marginal R2 (currentl nearly 43%) - need to redo without outlier

VCV_ES4 <- make_VCV_matrix(data = dat4, V = "lnRRV_ES", cluster = "Study_ID", obs = "ES_ID")
mod_ES6 <-rma.mv(yi = lnRR_ESa, V = VCV_ES4, mod = ~Stress_duration-1, random = list(~1|Study_ID,
                                                                                      ~1|ES_ID,
                                                                                      ~1|Strain),
                 test = "t",
                 data = dat4)
summary(mod_ES6) 
## 
## Multivariate Meta-Analysis Model (k = 87; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -41.2407   82.4815   92.4815  104.6947   93.2410   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0213  0.1459     27     no  Study_ID 
## sigma^2.2  0.0304  0.1744     87     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 85) = 293.4047, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 85) = 5.6275, p-val = 0.0051
## 
## Model Results:
## 
##                         estimate      se     tval  df    pval    ci.lb   ci.ub 
## Stress_durationAcute     -0.1285  0.1092  -1.1760  85  0.2429  -0.3457  0.0887 
## Stress_durationChronic    0.1644  0.0523   3.1420  85  0.0023   0.0603  0.2684 
##  
## Stress_durationAcute 
## Stress_durationChronic  ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES6) 
##   R2_marginal R2_coditional 
##     0.2226789     1.0000000
# Orchard plot 
orchard_plot(mod_ES6, mod = "Stress_duration", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Social enrichment

Does EE also include a manipulation of social environment (i.e., number of individuals in EE relative to control)?

VCV_ES5 <- make_VCV_matrix(data = dat5, V = "lnRRV_ES", cluster = "Study_ID", obs = "ES_ID")
mod_ES7<- rma.mv(yi = lnRR_ESa, V = VCV_ES5, mod = ~EE_social-1, random = list(~1|Study_ID,
                                                                                ~1|ES_ID,
                                                                                ~1|Strain),
                 test = "t",
                 data = dat5)

summary(mod_ES7)
## 
## Multivariate Meta-Analysis Model (k = 90; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -44.6612   89.3224   99.3224  111.7091  100.0541   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0237  0.1540     29     no  Study_ID 
## sigma^2.2  0.0283  0.1683     90     no     ES_ID 
## sigma^2.3  0.0385  0.1961      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 88) = 319.5169, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 88) = 1.0646, p-val = 0.3493
## 
## Model Results:
## 
##                      estimate      se    tval  df    pval    ci.lb   ci.ub 
## EE_socialNon-social    0.1083  0.1226  0.8838  88  0.3792  -0.1353  0.3519    
## EE_socialSocial        0.1791  0.1235  1.4501  88  0.1506  -0.0664  0.4246    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES7) 
##   R2_marginal R2_coditional 
##    0.01313274    0.58063226
# Orchard plot 
orchard_plot(mod_ES7, mod = "EE_social", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Exercise enrichment

Does the form of enrichment include exercise through a running wheel or treadmill?

mod_ES8<- rma.mv(yi = lnRR_ESa, V = VCV_ES, mod = ~EE_exercise-1, random = list(~1|Study_ID, 
    ~1|ES_ID,
    ~1|Strain),
     test = "t",
     data = dat)

summary(mod_ES8)
## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -45.5487   91.0975  101.0975  113.5965  101.8118   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0337  0.1837     30     no  Study_ID 
## sigma^2.2  0.0273  0.1653     92     no     ES_ID 
## sigma^2.3  0.0049  0.0702      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 90) = 309.5719, p-val < .0001
## 
## Test of Moderators (coefficients 1:2):
## F(df1 = 2, df2 = 90) = 1.9734, p-val = 0.1450
## 
## Model Results:
## 
##                         estimate      se    tval  df    pval    ci.lb   ci.ub 
## EE_exerciseExercise       0.1136  0.0762  1.4910  90  0.1395  -0.0378  0.2649 
## EE_exerciseNo exercise    0.1670  0.1002  1.6672  90  0.0989  -0.0320  0.3660 
##  
## EE_exerciseExercise 
## EE_exerciseNo exercise  . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES8) 
##   R2_marginal R2_coditional 
##   0.009818664   0.926129248
# Orchard plot 
orchard_plot(mod_ES8, mod = "EE_exercise", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

Age of enrichment

The age at which the individuals were exposed to environmental enrichment.
This needs to be redone after recategorising

mod_ES9 <- rma.mv(yi = lnRR_ESa, V = lnRRV_ES, mod = ~Age_EE_exposure-1, random = list(~1|Study_ID,
                                                                                       ~1|ES_ID,
                                                                                       ~1|Strain),
                  test = "t",
                  data = dat)

summary(mod_ES9) 
## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -42.9629   85.9259   97.9259  112.8577   98.9503   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0547  0.2340     30     no  Study_ID 
## sigma^2.2  0.0170  0.1305     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Residual Heterogeneity:
## QE(df = 89) = 299.7681, p-val < .0001
## 
## Test of Moderators (coefficients 1:3):
## F(df1 = 3, df2 = 89) = 3.0451, p-val = 0.0329
## 
## Model Results:
## 
##                          estimate      se     tval  df    pval    ci.lb   ci.ub 
## Age_EE_exposureAdult       0.1576  0.1030   1.5299  89  0.1296  -0.0471  0.3623 
## Age_EE_exposureJuvenile   -0.0146  0.1733  -0.0840  89  0.9332  -0.3589  0.3298 
## Age_EE_exposureUnclear     0.1693  0.0650   2.6053  89  0.0108   0.0402  0.2985 
##  
## Age_EE_exposureAdult 
## Age_EE_exposureJuvenile 
## Age_EE_exposureUnclear   * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
r2_ml(mod_ES9) 
##   R2_marginal R2_coditional 
##    0.03931516    0.99999999
# Orchard plot 
orchard_plot(mod_ES9, mod = "Age_EE_exposure", xlab = "lnRR", alpha=0.4) + 
  geom_errorbarh(aes(xmin = lowerPR, xmax = upperPR), height = 0, show.legend = FALSE, size = 1.1, alpha = 0.5) + # prediction intervals
  geom_errorbarh(aes(xmin = lowerCL, xmax = upperCL), height = 0.05, show.legend = FALSE, size = 2) + # confidence intervals
  geom_point(aes(fill = name),  size = 5, shape = 21)+ # mean estimate
  scale_size_continuous(range = c(1, 7))+ # change point scaling
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1.3), # border around the plot
        text = element_text(size = 24), # change font sizes
        legend.title = element_text(size = 15),
        legend.text = element_text(size = 13)) 

‘Pure effect sizes’

Comparing ES relative to stress we can see that the increase is about 18% but ES relative to control is 12% (i.e., this is the interaction).

We can also see that ES relative to control 9% so these individuals actually do better than baseline individuals. There is also no significant difference between ES and EC individuals (although estimate is slightly lower), meaning that enrichment in the presence of stress can result in learning and memory that is just as good as enriched individuals without stress

Enrichment relative to control

mod_E20 <- rma.mv(yi = lnRR_E2a, V = lnRRV_E2, random = list(~1|Study_ID, 
                                                          # ~ 1|Strain, does not run as we have NA
                                                          ~1|ES_ID),
                 test = "t",
                 data = dat)
summary(mod_E20) 
## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -84.0627  168.1253  174.1253  181.6579  174.4012   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0119  0.1089     30     no  Study_ID 
## sigma^2.2  0.0000  0.0000     92     no     ES_ID 
## 
## Test for Heterogeneity:
## Q(df = 91) = 18.3889, p-val = 1.0000
## 
## Model Results:
## 
## estimate      se    tval  df    pval   ci.lb   ci.ub 
##   0.1272  0.0619  2.0546  91  0.0428  0.0042  0.2501  * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
orchard_plot(mod_E20, mod = "Int", xlab = "lnRR")

Stress relative to control

mod_S20 <- rma.mv(yi = lnRR_S2a, V = lnRRV_S2, random = list(~1|Study_ID, 
                                                          # ~ 1|Strain, does not run as we have NA
                                                          ~1|ES_ID),
                 test = "t",
                 data = dat)
summary(mod_S20) #
## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -90.1569  180.3138  186.3138  193.8464  186.5897   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0167  0.1291     30     no  Study_ID 
## sigma^2.2  0.0000  0.0000     92     no     ES_ID 
## 
## Test for Heterogeneity:
## Q(df = 91) = 29.0591, p-val = 1.0000
## 
## Model Results:
## 
## estimate      se     tval  df    pval    ci.lb   ci.ub 
##  -0.1112  0.0660  -1.6859  91  0.0952  -0.2422  0.0198  . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
orchard_plot(mod_S20, mod = "Int", xlab = "lnRR")

Enrichment + stress relative to control

mod_ES20 <- rma.mv(yi = lnRR_ES2a, V = lnRRV_ES2, random = list(~1|Study_ID, 
                                                          # ~ 1|Strain, does not run as we have NA
                                                          ~1|ES_ID),
                 test = "t",
                 data = dat)
summary(mod_ES20) #
## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -88.6424  177.2848  183.2848  190.8174  183.5607   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0122  0.1103     30     no  Study_ID 
## sigma^2.2  0.0000  0.0000     92     no     ES_ID 
## 
## Test for Heterogeneity:
## Q(df = 91) = 22.5439, p-val = 1.0000
## 
## Model Results:
## 
## estimate      se    tval  df    pval    ci.lb   ci.ub 
##   0.0934  0.0639  1.4621  91  0.1472  -0.0335  0.2204    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
orchard_plot(mod_ES20, mod = "Int", xlab = "lnRR")

Enrichment + stress relative to stress

mod_E30 <- rma.mv(yi = lnRR_E3a, V = lnRRV_E3, random = list(~1|Study_ID, 
                                                               # ~ 1|Strain, does not run as we have NA
                                                               ~1|ES_ID),
                  test = "t",
                  data = dat)
summary(mod_E30) 
## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -91.7018  183.4037  189.4037  196.9362  189.6795   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0115  0.1072     30     no  Study_ID 
## sigma^2.2  0.0000  0.0000     92     no     ES_ID 
## 
## Test for Heterogeneity:
## Q(df = 91) = 26.3880, p-val = 1.0000
## 
## Model Results:
## 
## estimate      se    tval  df    pval   ci.lb   ci.ub 
##   0.1843  0.0629  2.9317  91  0.0043  0.0594  0.3092  ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
orchard_plot(mod_E30, mod = "Int", xlab = "lnRR")

Enrichment + stress relative to enrichment

mod_S30 <- rma.mv(yi = lnRR_S3a, V = lnRRV_S3, random = list(~1|Study_ID, 
                                                             # ~ 1|Strain, does not run as we have NA
                                                             ~1|ES_ID),
                  test = "t",
                  data = dat)
summary(mod_S30) 
## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##   logLik  Deviance       AIC       BIC      AICc 
## -82.4100  164.8200  170.8200  178.3526  171.0958   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.0000  0.0000     30     no  Study_ID 
## sigma^2.2  0.0000  0.0000     92     no     ES_ID 
## 
## Test for Heterogeneity:
## Q(df = 91) = 10.9070, p-val = 1.0000
## 
## Model Results:
## 
## estimate      se     tval  df    pval    ci.lb   ci.ub 
##  -0.0093  0.0503  -0.1855  91  0.8533  -0.1092  0.0906    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
orchard_plot(mod_S30, mod = "Int", xlab = "lnRR")

Combined orchard plot

#TODO need to figure out how to label y-axis. Also need to think about the most logical order to present these so that readers will understand what they are showing

mod_list <- list(mod_E20, mod_S20, mod_ES20, mod_E30, mod_S30)

mod_res <- lapply(mod_list, function(x) mod_results(x, mod = "Int"))

merged <- submerge(mod_res[[1]], mod_res[[2]],  mod_res[[3]], mod_res[[4]],  mod_res[[5]], mix = T)


#merged$mod_table$name2 <- c("EC vs CC", "CS vs CC", "ES vs CC", "ES vs EC", "ES vs CS")

orchard_plot(merged, mod = "Int", xlab = "lnRR", angle = 0) #+ geom_point(data = merged$mod_table, aes(y = name2, x = estimate))

Modelling with SMD

Stress

Intercept model

Why does the funnel plot look so strange and why is I2 higher?

mod_S0a <- rma.mv(yi = SMD_Sa, V = VCV_Sa, random = list(~1|Study_ID,
                                                         ~1|ES_ID,
                                                         ~1|Strain),
                test = "t",
                data = dat)
summary(mod_S0a) 
## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##    logLik   Deviance        AIC        BIC       AICc 
## -118.9545   237.9091   245.9091   255.9525   246.3742   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.6193  0.7869     30     no  Study_ID 
## sigma^2.2  0.4872  0.6980     92     no     ES_ID 
## sigma^2.3  0.0000  0.0001      6     no    Strain 
## 
## Test for Heterogeneity:
## Q(df = 91) = 27513.4051, p-val < .0001
## 
## Model Results:
## 
## estimate      se     tval  df    pval    ci.lb    ci.ub 
##  -0.3674  0.1666  -2.2044  91  0.0300  -0.6984  -0.0363  * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_S0a) 
##     I2_total  I2_Study_ID     I2_ES_ID    I2_Strain 
## 9.979242e-01 5.585262e-01 4.393980e-01 2.491240e-09
funnel(mod_S0a) 

funnel(mod_S0a, yaxis="seinv",) 

Environmental Enrichment

Intercept model

Why does the funnel plot look so strange and why is I2 higher?

mod_E0a <- rma.mv(yi = SMD_Ea, V = VCV_Ea, random = list(~1|Study_ID,
                                                         ~1|ES_ID, 
                                                         ~1|Strain),
                 test = "t",
                 data = dat)
summary(mod_E0a)
## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##    logLik   Deviance        AIC        BIC       AICc 
## -116.8436   233.6871   241.6871   251.7305   242.1522   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.3557  0.5964     30     no  Study_ID 
## sigma^2.2  0.5325  0.7297     92     no     ES_ID 
## sigma^2.3  0.0000  0.0000      6     no    Strain 
## 
## Test for Heterogeneity:
## Q(df = 91) = 14139.3073, p-val < .0001
## 
## Model Results:
## 
## estimate      se    tval  df    pval   ci.lb   ci.ub 
##   0.7058  0.1390  5.0791  91  <.0001  0.4298  0.9818  *** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_E0a)
##     I2_total  I2_Study_ID     I2_ES_ID    I2_Strain 
## 9.974155e-01 3.994255e-01 5.979900e-01 1.844209e-09
funnel(mod_E0a) 

Interaction

Intercept model

Why does the funnel plot look so strange and why is I2 higher?

mod_ES0a <- rma.mv(yi = SMD_ESa, V = VCV_ESa, random = list(~1|Study_ID,
                                                            ~1|ES_ID,
                                                            ~1|Strain),
                  test = "t",
                  data = dat)
summary(mod_ES0a)
## 
## Multivariate Meta-Analysis Model (k = 92; method: REML)
## 
##    logLik   Deviance        AIC        BIC       AICc 
## -128.8823   257.7646   265.7646   275.8080   266.2297   
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed    factor 
## sigma^2.1  0.8719  0.9338     30     no  Study_ID 
## sigma^2.2  0.5852  0.7650     92     no     ES_ID 
## sigma^2.3  0.0000  0.0001      6     no    Strain 
## 
## Test for Heterogeneity:
## Q(df = 91) = 12257.6601, p-val < .0001
## 
## Model Results:
## 
## estimate      se    tval  df    pval   ci.lb   ci.ub 
##   0.6295  0.1952  3.2255  91  0.0017  0.2418  1.0172  ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
i2_ml(mod_ES0a) #I2 is very high! 99.37%
##     I2_total  I2_Study_ID     I2_ES_ID    I2_Strain 
## 9.937219e-01 5.946336e-01 3.990883e-01 6.436368e-09
funnel(mod_ES0a)